############################################################################# 
# Prepared by Vladimir Chlouba 
# University of Richmond
# vchlouba@richmond.edu
# September 2025
#############################################################################
# Seeing is Believing: Voluntary Gender Quotas Change Female Leadership
# Stereotypes
# Replication R Script
#############################################################################

# loading the requisite R packages
library(ggplot2)
library(foreign)
library(data.table)
library(stargazer)
library(lme4)
library(arm)
library(memisc)
library(vegan)
library(MASS)
library(nnet)
library(spikeSlabGAM)
library(speedglm)
library(biglm)
library(ordinal)
library(erer)
library(mlogit)
library(scales)
library(nnet)
library(haven)
library(lfe)
library(sp)
library(raster)
library(lattice)
library(spdep)
library(cshapes)
library(countrycode)
library(margins)
library(ggmap)
library(dplyr)
library(lmtest)
library(sandwich)
library(fixest)
library(ivmodel)
library(dotwhisker)
library(interplot)
library(gridExtra)
library(effects)
library(ggtext) 


rm(list = ls())
#~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~

# loading the requisite dataset
namibia_constituency <- read.csv("chlouba_2025_R&P_replication_data.csv")

#~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~

namibia_constituency$event_time <- case_when(
  namibia_constituency$round == 3 ~ 0,   
  namibia_constituency$round == 5 ~ 1,   
  namibia_constituency$round == 6 ~ 2,   
  namibia_constituency$round == 7 ~ 3    
)

# both genders
mod.1 <- feols(mean_residual_female_leaders ~ i(event_time, SWAPO_average, ref = "0") | 
                 constituency + round, 
               data = namibia_constituency, 
               cluster = "constituency")

# men
mod.men <- feols(mean_residual_female_leaders_men ~ i(event_time, SWAPO_average, ref = "0") | 
                     constituency + round, 
                   data = namibia_constituency, 
                   cluster = "constituency")

# women
mod.women <- feols(mean_residual_female_leaders_women ~ i(event_time, SWAPO_average, ref = "0") | 
                       constituency + round, 
                     data = namibia_constituency, 
                     cluster = "constituency")

#~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~

# Figure 2

plot_data <- data.frame(
  Model = rep(c("both genders", "men", "women"), each = 3),
  Event_Time = rep(c(1, 2, 3), 3),
  Coefficient = c(
    mod.1$coefficients[1], mod.1$coefficients[2], mod.1$coefficients[3],  
    mod.men$coefficients[1], mod.men$coefficients[2], mod.men$coefficients[3],
    mod.women$coefficients[1], mod.women$coefficients[2], mod.women$coefficients[3]  
  ),
  SE = c(
    mod.1$se[1], mod.1$se[2], mod.1$se[3],
    mod.men$se[1], mod.men$se[2], mod.men$se[3],
    mod.women$se[1], mod.women$se[2], mod.women$se[3]
  )
)

# define jitter offsets to prevent overlapping whiskers
jitter_offsets <- c("both genders" = -0.1, "men" = 0, "women" = 0.1)

# define marker shapes and line types
marker_shapes <- c("both genders" = 16, "men" = 15, "women" = 17)  # circle, square, triangle
line_types <- c("both genders" = "solid", "men" = "dotdash", "women" = "dashed")

# create the plot
ggplot(plot_data, aes(x = Event_Time + jitter_offsets[Model], y = Coefficient, 
                      color = Model, shape = Model, linetype = Model)) +
  geom_point(size = 4) +
  geom_segment(aes(x = Event_Time + jitter_offsets[Model], 
                   xend = Event_Time + jitter_offsets[Model], 
                   y = Coefficient - 1.96 * SE, 
                   yend = Coefficient + 1.96 * SE), 
               linetype = "solid") +
  geom_line(aes(group = Model), size = 1) +
  scale_shape_manual(values = marker_shapes) +
  scale_linetype_manual(values = line_types) +
  theme_classic() +
  labs(title = " ",
       x = " ",
       y = "effect on support for female leaders") +
  
  # custom tick marks for survey years and political events
  scale_x_continuous(
    breaks = c(1, 1.5, 2, 2.3, 2.6, 3),
    labels = c(
      "<b>2012 survey</b>",
      "quota adopted<br>(June 2013)",
      "<b>2014 survey</b>",
      "elections<br>(Nov. 2014)",
      "new MPs enter<br>(March 2015)",
      "<b>2017 survey</b>"
    ),
    limits = c(0.8, 3.2)
  ) +
  
  # reference line at 0
  geom_hline(yintercept = 0, linetype = "dotted", color = "red") +
  
  # theme adjustments
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.text = element_text(size = 12),
    axis.text.x = element_markdown(size = 9, lineheight = 0.9),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12)
  )

#~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~
